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Over  the  past  six  months,  this  contract  has  funded  two  projects  in  full  and  one  project  in 
part.  The  two  fully  funded  projects  focus  on  the  application  of  random  coefficient  models  to 
wideband  high-resolution  direction  finding  and  transient  signal  detection  and  estimation.  The 
partially  funded  project  involves  the  analysis  of  nonlinear,  possibly  chaotic,  dynamical  systems. 
Progress  in  each  of  these  areas  is  described  below. 


Random  Coefficient  Models;  The  application  of  random  coefficient  models  to  narrowband 
high  ."esolution  direction  finding  has  been  very  successful  and  has  already  resulted  in  a  conference 
paper  to  be  presented  at  ICASSP92  (see  the  attached  paper  by  Jost  L  Williams).  It  has  been 
shown  that  the  random  coefficient  model  is  much  better  suited  to  modeling  sensor  array  data 
than  the  autoregressive  model  is.  A  simple  method  for  estimating  the  parameters  of  the  random 
coefficient  model  has  been  developed  and  applied  to  simulated  data.  Finally,  a  beamformer  for 
the  random  coefficient  model  has  been  developed  which  has  significantly  better  performance 
then  earlier  linear  predictive  beamformers.  A  journal  paper  describing  these  results  is  currently 
being  written  and  should  be  submitted  to  the  IEEE  Transactions  on  Signal  Processing  by  this 
Spring.  There  is  also  progress  being  made  toward  the  final  goal  of  applying  these  techniques  to 
wideband  direction  finding. 


Transient  Signal  Detection:  The  application  of  the  wavelet  transform  to  the  detection 
of  transient  signals  with  an  array  of  sensors  is  being  examined.  This  approach  has  led  to  a 
directional  multirate  filter  bank  structure  that  decomposes  the  incoming  signal  into  decaying 
exponentials.  This  filter  bank  is  also  capable  of  adapting  towards  an  improved  estimate  of  the 
structure  of  the  transient  signal  which,  consequently,  also  improves  the  detection  performance. 
Publication  of  these  results  will  proceed  after  simulations  and  comparisons  to  other  transient 
signal  detectors  are  complete. 

Nonlinear  System  Identification:  System  identification  algorithms  that  depend  on  gradient 
descent  methods  have  been  found  to  degrade  significantly  if  the  time-series  or,  equivalently,  the 
system  that  produced  the  time-series  is  chaotic  (see  the  attached  paper  by  Drake  k.  Williams). 
A  careful  analysis  of  these  degradations  has  led  to  algorithms  which  are  much  less  sensitive  to 
the  potentially  chaotic  nature  of  these  nonlineau’  systems.  Analysis  of  these  systems  and  their 
time-series  continues  with  the  eventual  goal  being  a  very  general  system  identification  algorithm 
to  be  applied  to  the  time-series  produced  by  nonlinear  discrete-time  systems. 
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To  be  presented  at  ICASSP,  San  Francisco,  California,  March  1992. 
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ABSTRACT 

This  paper  applies  the  random  coefficient  model  to  array 
processing,  specificaUy  in  the  design  of  a  beamformer  for  di¬ 
rection  finding.  This  model  is  similar  to  the  autoregressive 
(AR)  model,  except  the  coefficients  are  allowed  to  change 
with  time  instead  of  remaining  constant;  thus  allowing  the 
beamformer  to  better  model  any  additive  noise  or  signal  cor¬ 
relations  in  the  observations.  Through  the  use  of  a  binary 
hypothesis  test,  it  is  shown  that  random  coefficient  models 
better  fit  typical  array  data  than  do  AR  models.  A  Kalman 
filter  is  presented  that  has  the  array  observations  as  inputs 
and  the  pairametets  of  the  ramdom  coefficient  model  as  out¬ 
puts.  A  new  beamformer  based  on  the  random  coefficient 
model  is  derived  that  is  similar  to  the  constant  coefficient 
linear  predictive  (CCLP)  beamformer.  The  two  beamform- 
ers  are  compared  and  it  is  shown  that  the  random  coefficient 
beamformer  outperforms  the  CCLP  beamformer. 


1.  INTRODUCTION 

The  linear  model  is  probably  the  most  popular  model  in  use 
in  engineering  and  science  because  it  is  simple  and  yet  pow¬ 
erful.  Many  times,  computational  advantages  due  to  the 
simplicity  of  the  linear  model  far  outweigh  any  performance 
gain  achieved  by  more  complicated  models.  The  commonly 
used  autoregressive  (AR)  model  is  one  such  linear  model. 
In  array  processing  the  AR  model  leads  to  linear  predictive 
(LP)  beamforming  [1]  where  the  output  of  a  selected  sensor, 
say  the  moth,  is  estimated  as  a  weighted  linear  combination 
of  the  other  sensor  outputs.  Assuming  narrowband  signals 

and  letting  Xm(t),t  =  0 . A  —  1,  be  snapshots  of  the 

Fourier  transform  of  the  mth  sensor’s  output  at  the  signals’ 
common  center  frequency,  the  LP  model  for  X  fFlQ  (t)  is 

Xmo(t)  =  -  51  (1) 


where  {om}  is  a  set  of  complex-valued  weights  to  be  found. 

When  compared  to  other  beamformers  with  roughly 
equivalent  computational  complexity,  such  as  the  minimum 
variance  distortionless  receiver,  the  LP  method  possesses 
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superior  resolution  properties  [1].  However,  the  perfor¬ 
mance  of  the  LP  beamformer  is  highly  dependent  on  both 
the  signals’  directions  of  propagation  and  the  selection  of 
mo.  Also,  in  many  cases  the  direction-of-arrival  estimation 
bias  is  quite  high.  The  primary  causes  of  this  estimation 
bias  are  additive  noise  in  the  observations  and  correlations 
between  propagating  signals.  Additive  noise  and  signal  cor¬ 
relations  do  not  obey  an  AR  relationship  across  the  array; 
so  the  model  in  (1)  is  no  longer  accurate. 

In  order  to  keep  the  simplicity  of  the  linear  model  of  (1) 
and  yet  more  accurately  match  the  array  data,  the  random 
coefficient  model  is  proposed.  This  model  is  identical  to 
the  AR  model  except  the  weights  Om  are  random  instead  of 
constant  and  (1)  becomes 

-  53  (2) 

mfimo 

The  coefficients  can  be  expressed  as 

nm(k)  =  firn  ^  r^m(k),  (3) 


where  fim  is  the  mean  of  the  coefficients  for  the  mth  sensor 
and  the  i'm(t)  are  zero-mean  independent  identically  dis¬ 
tributed  random  variables.  The  random  coefficient  model 
of  (2)  keeps  the  simple  linear  form  of  the  AR  model  and, 
hence,  its  computational  advantages.  Since  the  constant 
coefficients  of  an  AR  model  can  be  expressed  in  terms  of 
(3)  with  the  I'm  (It)  equal  to  zero,  the  random  coefficient 
model  can  actuaUy  be  thought  of  as  a  generalization  of  the 
AR  model.  Therefore,  the  random  coefficient  model  will  fit 
the  array  data  at  least  as  well  as  the  AR  model  and  will 
simplify  to  the  AR  model  when  it  is  the  better  choice. 


2.  TESTING  OF  THE  MODELS 

Through  the  use  of  a  binary  hypothesis  test  developed  by 
Breusch  and  Pagan  [2],  it  is  possible  to  test  which  model, 
AR  or  random  coefficient,  will  better  fit  typical  array  data. 
The  test  indicates  that  if  the  residuals  from  an  AR  model 
have  a  constant  variance,  an  AR  model  is  better  suited; 
otherwise,  a  random  coefficient  model  should  be  considered. 
The  test  is  performed  by  generating  a  Lagrange  multiplier 
(LM)  statistic  from  generated  array  data.  If  the  nuU  hy¬ 
pothesis  of  an  AR  model  is  true,  then  the  LM  statistics 
will  be  asymptotically  distributed  as  chi-squared  with  M 
degrees  of  freedom,  where  M  is  the  number  of  sensors  in 
the  array.  The  methpd  used  for  this  paper  was  to  generate 
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Figure  1:  Reaulu  of  the  hypotheaii  test.  The  ten  degree  of 
freedom  chi-squered  distribution  is  the  solid  line  while  the  nor¬ 
malised  histogram  of  the  LM  statistics  is  the  broken  line.  The 
array  data  was  from  a  uniform  linear  array  of  10  sensors  with 
spacing  -j.  There  are  two  signals  present  at  -10*  and  30*,  the 
noise  is  additive  white  Gaussian  and  the  SNR  is  6  d£.  There  are 
100  independent  samples  being  used  {N  —  100). 

array  data  for  a  given  number  of  signals  at  their  respective 
ang^les  of  arrival,  and  a  given  number  of  sensors,  time  sam¬ 
ples,  and  signal-to-noise  ratio.  The  LM  statistic  given  in  [2] 
was  calculated  for  a  large  number  of  sets  of  data  and  the 
normalized  histogram  of  these  LM  statistics  was  calculated 
and  plotted  against  a  chi-squared  density  function  with  M 
degrees  of  freedom  to  see  how  closely  the  two  matched.  Fig¬ 
ure  1  shows  an  example  of  this  test  in  which  the  normalized 
histogram  of  the  LM  statistics  clearly  does  not  match  the 
chi-squared  density  plot.  In  fact,  it  was  found  that  the 
AR  hypothesis  is  rejected  for  even  the  simplest  cases  if  any 
additive  noise  or  signal  correlation  is  present.  Therefore, 
it  would  seem  that  the  random  coefficient  model  is  better 
suited  for  array  processing  applications. 

3.  RANDOM  COEFFICIENT  BEAMFORMER 

As  the  random  coefficient  model  has  been  shown  to  fit  the 
data  more  accurately  than  the  AR  model,  a  beamforming 
algorithm  based  on  this  model  has  been  developed.  The 
remainder  of  the  paper  will  present  the  algorithm  for  esti¬ 
mating  the  random  coefficients  and  then  the  derivation  of 
the  random  coefficient  beamformer. 

3.1.  Estimating  the  Coefficients 
The  first  step  is  to  estimate  the  random  coefficients.  This 
is  performed  with  the  Kalman  filter  outlined  in  [3]  and  [4]. 
The  following  set  of  vectors  is  introduced  to  ease  the  com¬ 
putations; 

a'(fc)  =  [oo(t)-o„._,(t)o„,+,(it)--ow_,(fc)] 

b'  =  fimo+i’-  Pu-i] 

y'(k)  =  [^o(h)  + 

z'{k)  =  [Xo(it)-A-„,-,(fc)Xm.+i(t)-A-M-i(h)] 

where  M  is  the  number  of  sensors  and  (')  denotes  Eermitian 
transpose.  An  error  covariance  matrix  is  defined  as 

^  =  k?>l.  «,;  =  0, . . . ,  mo  -  1,  mo  -f-l, . . . ,  Af  -  1 

where  (1)}  =  6{k  —  Ftom  (3)  and  [4]  the  es¬ 

timates  of  a(it)  are  calculated  from  the  K^man  filter  equa¬ 
tions 


£(i/i  -  1)  =  E 

k{k/k)  =  b  -I-  E(i/i  -  l)z{k)h-\k) 

S(k/k)  =  L(k/k-l)- 

t(k/k  -  l)z{k)h-^(k)z'(k)t{k/k  -  1)  (4) 

where  k(k/k)  is  the  estimate  of  a(h)  at  time  k  given  all  the 
observations  up  to  k,  tl(k/k)  is  the  error  covariance  matrix 
estimate  of  a(h)  at  time  k  given  all  the  observations  up  to 
k,  and  is  the  noise  variance  estimate. 

To  implement  the  Kalman  filter,  initial  estimates  of  b,  E, 
and  cr’  must  be  calculated.  The  authors  will  present  a 
method  that  calculates  the  least-squares  estimates  of  these 
terms  following  [4].  Introducing  a  noise  term  into  (2)  yields 

A-„.{t)  =  -  ^  a„(ib)A-„(fc)H-n„o(i:)  (5) 

wnere  nmo  is  white  Gaussian  noise  and  is  uncorrelated 
with  the  observations  Xm{k).  The  noise  associated  with 
an  arbitrary  sensor,  say  the  mth,  has  the  property 
E{nm(»nm(t)}  =  ^0  -  Using  (3),  equation  (5)  can 
be  written  as 

Xmo(ifc)  =  - 

mf<mo  trip 

+  nmo(*;) 

=  -Y  +  (6) 

where 

«(fc)  =  nmo(h)-  Y^  t'm{t)A’m(i).  (7) 

fnp 

Defining 

u'  =  [tt(0)  ■  •  •  «(A' —  1)] 

y'  =  [X„,{G)-X,r.,{N-\)] 

Z'  =  [z(0)...z(fV-l)], 

equation  (6)  can  be  written  as 

y  =  — Zb-i-  u. 

Using  ordinary  least  squares,  the  estimate  of  b  is 

b= -(Z'Z)"*Z'y.  (8) 

Taking  the  expected  value  of  the  magnitude  squared  of 
(7)  yields 

r{|«(t)|»}  =  £{|  «„.(t)  -  Y  I*} 

which  reduces  to 


i(k/k  -  1)  =  b 


(9) 


where  the  other  terms  &re  zero  because  the  noise  is  nncorre- 
Uted  with  the  signals  and  coefficients.  Assuming  the  Xi(k) 
and  ai(k)  are  zero  mean  and  Gaussian,  (9)  can  be  expressed 

as 

•^"»o  i¥mo 

£{u.(k)x.(t)}£{i.;(t)x;(t)}  + 

£{^.(k)x;(i)}£{,.;(k)x:(k)}] 

which  reduces  to 


The  right-hand  expected  value  can  be  expanded  as  in  (9) 

=  ^^£{a^a;}£{X.Xj%,}  +  A(/)  +  B(l)  (14) 
•  j 

where 

i  J 


£{|  u(k)  |»}  =  (r*+  53  52  £{‘'.(t)>'j(imx^(k)x;(k)}. 

•flmo  Jfimo 

(10) 

Since 

i«»(t  +  oi*=r{i«(t+on  +  e<  (11) 

where  the  is  a  zero  mean  error  with  unknown  variance, 
(10)  can  be  used  in  (11)  as 

|«(i  +  /)|*  =  ^»  + 

52  + o^;(i + 1))+(< 

ill  mg  jglmg 

or  in  vector  notation  as 

e  a  Xs  4  r 


for  {  =  — (Af  — 1), . .  The  expected  values  in  (14)  can 

be  estimated  by  averaging  over  the  available  data  and  co¬ 
efficient  estimates,  e.g.,  £{a’Xi}  =  ^  kj(k)Xi(k) 

The  spatial  Fourier  transform  of  (14)  can  be  calculated 
for  any  array  geometry,  but  if  a  linear  uniform  array  is 
used,  the  transform  can  be  parameterized  in  terms  of  just 
the  direction  of  look  9.  Therefore,  the  random  coefficient 
beamformer  for  a  uniform  linear  array  steered  in  direction 
9  can  be  written  as 


<r»  ~>1(9)-B(9) 


M  M 


(15) 


5252£;{o.«;}e-^»'^^- 


imi  jml 


where  e'  *=  (|tt(0)|*  •  ■  •  lu(Ar  —  l)|*]ir*  =  =  where  d  is  the  spacing  between  sensors,  A  is  the  wavelength 

^i.j  ■  •  •  <^2,1  •  •  •^Ar-i.v-il>  *^<1  lit*  propagatmg  signals,  and  A  and  B  are  the  spatial 


X  = 


1  x,(o)x;(o)  x,(o)x;(o)  •••  x,(o)xj;,_,(o)  Xs(o)xr(o) 

1  XiiN-l)X;(N-l) 


XM~i{o)x-^_,{o) 
Xu.r{N-l)X’u.,{N-l)  . 


Using  ordinary  least  squares,  the  estimate  of  s  is 

8  =  (X'X)-‘X'e,  (12) 

bom  which  the  estimates  of  and  £  are  extracted. 

3.2.  The  Beamformer 

Once  the  random  coefficients  have  been  estimated,  they  can 
be  used  in  a  beamformer  to  calculate  the  direction(8)  of  the 
measured  8ignal(8).  Following  the  derivation  of  the  constant 
coefficient  LP  beamformer,  (5)  can  be  rewritten  as 

nmo(h)  =  Xmo(i)4  5^  Om(k)Xm(k) 

A/-1 

=  52  am(*)Xm(k)  (13) 

msO 

if  8mo(k)  =  1  for  t  =  0, ...,iV  —  1.  Multiplying  (13)  with 
a  similar  expression  for  n^^^,(h)  and  taking  the  expected 
value  gives  (dropping  the  time  index  for  brevity) 

=  f(iyn  =  5252^(»‘*rA’lX;+,}. 


Fourier  trajisforms  of  j4(/)  and  B(l),  respectively,  e.g., 

JV-3 

A(0)=  51  A(l)e-'*’'<'""*. 

I— («-i) 

3.3.  Simulation  Results 

The  random  coefficient  beamformer  of  (15)  can  be  com¬ 
pared  to  a  constant  coefficient  LP  (CCLP)  beamformer. 
A  typical  example  of  this  comparison  is  shown  in  figure  2 
where  the  predictive  element  mo  is  the  first  sensor  for  both 
beamformers.  It  is  seen  that  while  both  have  sharp  peaks 
at  —10*  and  30*,  the  random  coefficient  beamformer  has  a 
lower  noise  floor  with  none  of  the  spurious  peaks  that  might 
be  confused  as  signals.  This  improvement  is  due  to  the  ran¬ 
dom  coefficient  model  being  able  to  fit  the  noisy  array  data 
better  than  the  AR  model.  As  was  expected  from  the  re¬ 
sults  of  the  hypothesis  test,  in  the  comparisons  run  by  the 
authors  the  random  coefficient  beamformer  clearly  outper¬ 
forms  the  CCLP  almost  every  time  and  always  performs  at 
least  as  well  as  the  CCLP. 


Figure  2:  Comparison  of  the  CCLP  and  random  coefficient  beamformers.  The  arra)’  data  used  was  the  same  as  that  described  in 
figure  1. 


4,  CONCLUSION 

The  random  coefficient  model  was  introduced  for  use  in 
array  processing.  It  was  shown  through  the  use  of  a  hy¬ 
pothesis  test  that  the  random  coefficient  model  will  better 
fit  typical  array  data  than  the  constant  coefficient  autore¬ 
gressive  model.  A  beamforraer  using  the  random  coefficient 
model  was  developed  (15)  as  well  as  the  method  for  esti¬ 
mating  the  random  coefficients  using  a  Kalman  filter.  The 
random  coefficient  beamformer  was  then  compared  to  the 
familiar  constant  coefficient  linear  predictive  beamformer 
and  was  shown  to  be  a  dramatic  improvement. 
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ABSTRACT 

The  extreme  Mnsitivity  of  a  chaotic  system’s  steady  state 
response  to  small  changes  in  its  initial  conditions  makes  long 
term  prediction  of  the  evolution  of  such  a  system  difficult,  if 
not  impossible.  In  the  framework  of  parameter  estimation, 
we  show  how  this  sensitivity  can  hinder  attempts  to  deter¬ 
mine  model  parameters  that  will  reproduce  a  target  chaotic 
time  sequence.  Specifically,  a  waveform  error  minimization 
technique  based  on  gradient  descent  optimization  is  not  weU 
suited  for  estimating  the  parameters  of  a  strongly  chaotic 
system.  We  propose  a  modification  of  this  minimization 
procedure  that  avoids  some  of  the  obstacles  present  when 
estimating  the  parameters  of  a  chaotic  system. 

1.  INTRODUCTION 

Chaos — unpredictable  deterministic  behavior — has  been  ob¬ 
served  in  phenomena  ranging  from  chemical  reactions  [l]  to 
solar  flares  [2].  Modelling  time  sequences  derived  from  such 
processes  can  provide  insight  into  the  underlying  physics 
that  drive  them.  Unfortunately,  the  intrinsic  sensitivity  of 
chaotic  systems  makes  them  difflcult  to  model;  a  represen¬ 
tation  with  enough  freedom  to  correctly  reproduce  chaotic 
behavior  will  itself  be  extremely  susceptible  to  small  varia¬ 
tions  in  its  parameters. 

The  realization  that  long-term  prediction  of  certain  com¬ 
pletely  deterministic  systems  was  impossible  sparked  inter¬ 
est  in  a  new  area  of  Dynamical  Systenu,  an  area  dealing 
with  the  phenomenon  of  chaos.  The  classic  description  of  a 
chaotic  system  usually  includes  the  phrase  “sensitive  depen¬ 
dence  on  initial  conditions”  [3].  Sensitive,  in  this  context, 
refers  to  the  exponential  rate  at  which  initially  close  tra¬ 
jectories  on  the  attractor  diverge.  This  sensitivity  can  be 
quantified  by  the  spectrum  of  Lyapunov  exponents  associ¬ 
ated  with  the  attractor. 

Quatieri  and  Hobtetter  [4]  wished  to  determine  the  pa¬ 
rameters  and  initial  conditions  of  a  nonlinear  difference 
equation  whose  solution  would  be  as  close  as  possible  to 
some  target  time  sequence  generated  by  a  dynamical  sys¬ 
tem.  They  derived  a  gradient  descent  method  that  mini¬ 
mized  the  waveform  error  between  the  solution  of  the  dif¬ 
ference  equation  model  and  the  target  sequence. 

We  will  show  that  the  waveform  error  surface  is  not  well- 
behaved  if  the  target  sequence  is  generated  by  a  chaotic  sys- 

Thit  work  was  supported  in  psirt  by  the  Office  of  Naval  Re¬ 
search  under  contract  N00014-91-J-4129. 


tern.  In  particular,  in  the  neighborhood  of  the  global  mini¬ 
mum  at  least  one  eigenvalue  of  the  Hessian  of  the  waveform 
error  increases  exponentially  as  a  function  of  the  length 
of  the  target  sequence.  This  unbounded  growth  sets  the 
global  minimum  at  the  bottom  of  a  deep  trench,  rendering 
gradient  descent  techniques  impractical.  Consequently,  we 
have  modified  the  waveform  error  minimization  procedure 
by  taking  into  account  the  behavior  of  the  error  surface  as  a 
function  of  target  sequence  length.  This  improved  optimiza¬ 
tion  technique  provides  a  much  wider  basin  of  attraction  for 
the  global  minimum  than  the  original  method. 

2.  DYNAMICAL  SYSTEMS  AND  CHAOS 

Suppose  h:MxR^-»Misa  parameterized  discrete- time 
dynamical  system  defined  on  a  smooth  compact  manifold 
M  such  that 

y[n]  =  h(y[n-l],p).  (1) 

We  assume  that  this  system  is  stable,  and  further  that  the 
state  y[n]  converges  onto  an  attractor  A  C  M  as  n  tends 
to  infinity  for  any  initial  condition  y[— 1]  contained  in  the 
basin  of  attraction  of  A. 

It  can  be  shown  that  the  variation  of  the  state  y[n]  with 
respect  to  initial  conditions  y[— 1]  is  given  by 

n— 1 

D»[-iiyW=  (2) 

where  the  operator  Di  applied  to  the  vector-valued  func¬ 

tion  f(x)  results  in  a  matrix  with  elements  (D^f)^  ^  = 
dfi{x)/dxf.  The  Lyapunov  exponents  quantify  the  average 
rate  at  which  small  perturbations  of  the  initial  condition 
are  exponentially  amplified  or  attenuated  upon  iterations 
of  the  system  [5].  An  infinitesimal  deviation  dy[-l]  will 
result  in  a  deviation 

n— 1 

dy[n]  =  JJ  D,h(y(i  -  l].p)dyl-l].  (3) 

&nd  for  almost  all  dy[— 1] 

IMyWlI  «  e*<"+»||dy[-l]|l.  (4) 

where  X  is  the  largest  Lyapunov  exponent  of  h  on  A.  Eq  (4) 

is  equivalent  to  saying  that  ^^»^(y[*  ~  l]iP) 


eigenvalue  that  grows  on  average  and  in  absolute  value  as 
A  system  is,  by  definition,  chaotic  if  it  has  least 
one  positive  Lyapunov  exponent,  indicating  its  exponential 
sensitivity  to  small  variations  in  initial  conditions. 

Similarly  one  can  show  that  the  dependence  of  the  state 
on  small  variations  of  the  system’s  parameters  is  given  by 


vectors  (x[n]),  where  x[n]  *  (z(f»],*[f>-l) . *[n-m+l]), 

with  the  embedding  dimension  m  snitably  chosen  .  The 
diifeomorphic  relationship  between  the  true  trajectory  in 
phase  space  and  the  reconstructed  one  preserves  certain 
quantities,  namely  the  Lyapunov  exponents.  Therefore  the 
dynamical  system 


n-l  /  n-l 

i>»-yW=  ^  (  n 

The  term  Dph  converts  small  deviations  in  the  parameters 
into  small  deviations  in  the  state  which  are  then  propagated 
forward  by  the  product  Dph.  This  coupling  between  pa¬ 
rameters  and  state  implies  that  a  chaotic  system  will  be 
extremely  sensitive  not  only  to  variations  in  its  initial  con¬ 
ditions  but  to  variations  in  its  parameters  as  well. 


j  Dph(y[i],p).  (5) 


3.  WAVEFORM  ERROR  MINIMIZATION 

Suppose  we  have  a  scalar  time  sequence  z[n]  =  v(y[n])  de¬ 
rived  from  a  dynamical  system  via  v  :  M  —»  R.  We  assume 
this  sequence  is  the  solution  of  an  m*'*  order  nonlinear  dif¬ 
ference  equation  with  a  known  form,  but  depending  on  k 
unknown  parameters  p.  An  estimate  p  of  these  parameters 
produces  the  time  sequence  estimate 

x(n]  =  /(x[n  —  1],  p)  with  0  <  n  <  .A^,  (6) 

where  x(t»  —  1)  =  (z[n  —  1],  i[n  —  2], ... ,  i[n  —  m])^  is  the 
vector  of  the  last  m  values  of  z  at  time  n.  For  simplicity  we 
assume  that  the  initial  conditions  x[— 1]  are  known  exactly' ; 
let  x[-l]  =  x(— 1].  We  wish  to  find  the  parameters  that 
minimize  the  waveform  error 

N-3 

= -^  ]^(z[n]-z[n])’.  (7) 

nvO 

Quatieri  and  Bofstetter  use  a  gradient  descent  method 
to  minimize  the  waveform  error  with  respect  to  parameters. 
An  initial  estimate  p  of  the  parameter  values  is  iteratively 
updated 

p  —  p-  fi^DpEn)^  (8) 

so  that  the  error  decreases  at  each  step.  The  step  size  ft 
is  chosen  so  that  the  error  decreases  at  each  iteration;  the 
minimization  procedure  terminates  when  the  waveform  er¬ 
ror  falls  below  a  specified  threshold. 

In  order  to  better  understand  the  behavior  of  the  error 
surface,  we  expand  En  about  the  global  minimum  p; 


x[«]  =  gWn-l],p)  = 


/(x[n-  l],p)  ■ 
x[n  -  1] 
x[n  —  2] 


x[f»  -  m  -f  1] 


(10) 


has  the  same  Lyapunov  exponents  as  the  original  dynami¬ 
cal  system.  The  gradient  of  the  scalar  time  sequence  (z[n]) 
is  simply  the  first  row  of  the  matrix  Dpx[n].  If  the  dy¬ 
namical  system  that  produced  the  sequence  is  chaotic,  then 
||Dpz[n]||  will  grow  exponentially  fast  with  increasing  n, 
since  it  generally  won’t  be  orthogonal  to  the  eigenvector 
along  which  the  exponential  expansion  is  taking  place.  Thus 
the  Hessian  of  the  waveform  error  in  Eq  (9),  composed  of 
a  sum  of  outer  products  of  the  vectors  (Dpz[n]),  has  an 
increasingly  large  eigenvalue.  As  N  increases  the  global 
minimum  will  become  sandwiched  between  two  increasingly 
steep  walls —  not  ideal  conditions  for  gradient  descent  op¬ 
timization. 

As  will  be  seen  in  the  next  section,  the  waveform  error 
seems  well-behaved  for  short  chaotic  sequences;  the  expo¬ 
nential  amplification  of  parameter  mismatch  has  little  time 
over  which  to  markedly  modify  the  sequence  estimate.  Our 
modification  of  the  waveform  error  minimization  procedure 
takes  advantage  of  this  phenomenon.  Instead  of  trying 
to  optimize  our  parameter  estimates  for  the  whole  target 
sequence  at  once,  we  sequentially  minimize  the  waveform 

errors  Eo,Ei . Es-  Once  En  sinks  below  some  fixed 

threshold,  we  repeat  the  minimization  process  on  En^i ,  us¬ 
ing  the  last  estimates  of  the  parameters  as  the  initial  guess 
for  the  next  step.  Since  each  error  surface  has  the  global 
minimum  in  common,  successive  minimization  of  the  errors 
forces  the  parameter  estimates  closer  to  their  true  value. 
We’ve  effectively  expanded  the  global  minimum’s  basin  of 
attraction  to  that  of  a  length-one  sequence,  independent  of 
the  true  sequence  length. 


4.  EXAMPLES 

A  simple  system  capable  of  exhibiting  chaotic  behavior  is 
the  logistic  equation 


£jv  »  -^dp^  ^  ^  (Dpi[n])’’  (Dpz[n])j  dp,  (9) 

where  dp  represents  an  infinitesimal  deviation  from  the  true 
parameters.  A  remarkable  result  by  Takens  [6]  states  that 
under  the  proper  conditions,  a  scalar  time  sequence  can  be 
‘time  delay  embedded’  into  R"*,  revealing  a  diffeomorphic 
copy  of  the  phase  space  dynamics  that  generated  the  se¬ 
quence.  The  embedding  is  represented  by  the  sequence  of 

'  The  origin  of  the  sequence  can  always  be  shifted  to  the  right 
by  m  samples. 


x[n]  =  p,x[n-l](l-x[n-l]),  (11) 

which  is  both  a  scalar  dynamical  system  and  a  first  order 
nonlinear  difference  equation.  The  model  exhibits  markedly 
different  types  of  steady  state  behavior  depending  on  the 
choice  of  Pi .  A  parameter  value  of  pi  =  3.5  results  in  a 
period-four  oscillation;  in  contrast,  a  value  pi  =  3.7  pro¬ 
duces  chaotic  behavior.  Figure  1  contrasts  the  waveform 
error  for  both  cases,  for  sequences  of  length  N  =  and 
an  initial  condition  z[— 1]  =  0.42.  While  relatively  flat 
and  smooth  in  the  periodic  case,  the  waveform  error  in  the 
chaotic  case  is  riddled  with  local  minima  and  the  global 


Figure  1;  The  waveform  errors  of  periodic  (above)  and 
chaotic  (below)  target  sequences  of  length  87,  generated  by 
the  logistic  equation. 


minimum  has  the  very  narrow  basin  of  attraction  as  dis¬ 
cussed  above. 

Figure  2  illustrates  the  behavior  of  the  waveform  er¬ 
ror  for  the  chaotic  logistic  system  with  respect  to  the  se¬ 
quence  length.  As  noted  previously,  error  surfaces  for  short 
sequence  lengths  are  relatively  smooth,  giving  parameter 
estimates  with  relatively  large  errors  a  greater  chance  of 
converging  to  the  true  parameter. 

Figure  3  compares  the  performance  of  the  originrJ  wave¬ 
form  error  minimization  technique,  which  tries  to  minimize 
the  waveform  error  En  directly,  and  our  modified  version. 
As  expected  from  the  general  appearance  of  the  waveform 
error,  an  initial  deviation  of  only  5  x  10~*  in  the  model 
parameter  gets  trapped  almost  immediately  in  a  local  min¬ 
imum,  and  the  error  never  descends  below  the  specified 
threshold  of  10~*.  Our  extension  method,  on  the  other 
hand,  correctly  identifies  the  true  parameter  after  extend¬ 
ing  the  target  sequence  to  N  =  87,  even  though  the  initial 
deviation  from  the  true  parameter  value  was  5  x  10~^;  six 
orders  of  magnitudes  larger.  In  fact,  any  initial  parame¬ 
ter  value  within  the  logistic  equation’s  usual  working  range 
po  €  [0, 4]  will  be  converge  to  the  the  true  parameter  value. 

An  example  of  a  two-parameter  chaotic  system  is  the 
Henon  system 

I»i[n]  =  1-P>»?[n-1]  +  P3[n-1]  (12) 

(13) 

with  parameter  values  pi  =  1.4,  ps  =  0.3.  If  we  consider 
the  time  sequence  produced  by  the  first  variable  (viN)  out 
difference  equation  model  has  the  form 

*[n]  =  1 -po**[n  -  1]  +  Pi*(n- 2],  (14) 

with  initial  conditions  z[— 1]  —  |/i[— 1]  =  0.948586  and 
*1-2]  =  wl-lj/pj  =  0.425317. 


Fignre  2:  The  waveform  error  as  a  function  of  target  se¬ 
quence  length.  In  this  case  the  sequence  was  generated  by 
ue  chaotic  lo^tic  equation.  The  error  surface  is  relatively 
smooth  and  snaUow  for  short  sequences,  and  becomes  in¬ 
creasingly  rough  as  more  of  the  sequence  is  considered. 


Figure  4  shows  that  our  extension  method  outperforms 
the  original  waveform  error  minimization  technique.  Initial 
parameter  estimates  with  errors  of  more  that  10~’  are  re¬ 
duced  by  six  orders  of  magnitude.  The  original  method, 
initiated  with  deviations  of  only  5  x  10~*  in  both  parame¬ 
ters,  is  immediately  trapped  in  a  local  minimum. 

Unlike  the  one  dimensional  case,  for  this  two-parameter 
system  our  method  does  not  produce  estimates  that  con¬ 
verge  to  the  true  parameter  values.  Figure  5  shows  the 
error  surface  £so  in  a  small  neighborhood  of  the  global 
minimum.  As  expected,  the  sharp  gradient  discussed  in 
previous  sections  is  in  evidence.  However,  there  also  seems 
to  be  a  continuous  range  of  parameter  values  that  generate 
the  same  waveform  as  those  located  at  the  global  minimum. 
This  alignment  is  representative  of  the  true  behavior  of  the 
Henon  system  and  is  not  an  artifact  of  the  conversion  from 
dynamical  system  to  difference  equation.  An  examination 
of  the  gradient  of  (yi[n])  with  respect  to  the  parameters 
shows  that  while  the  (DpVi  [n])  grow  quickly  for  increasing 
n  as  expected,  they  also  tend  to  align  themselves  along  a 
common  axis.  The  sum  of  outer  products  in  Eq.  9  is  domi¬ 
nated  by  the  matrices  formed  from  these  increasingly  large 
gradients  that  all  point  in  the  same  direction,  and  conse¬ 
quently  the  Hessian  appears  singular.  This  behavior  is  not 
typical  of  chaotic  dynamical  systems  in  general. 

6.  CONCLUSION 

Using  concepts  from  the  discipline  of  Dynamical  Systems 
we  have  shown  how  the  sensitive  dependence  of  a  chaotic 
system  on  its  initial  conditions  can  induce  an  analogous  de¬ 
pendence  on  its  parameters.  Takens’  embedding  theorem 
allowed  us  to  transplant  the  phase  space  based  notion  of 
Lyapunov  exponents  which  quantify  this  sensitive  depen- 


Figure  3:  Comparison  of  waveform  error  minimization 
tedniques  using  the  chaotic  logistic  equation.  Even  good 
initial  parameter  estimates  get  trapped  in  local  minima 
when  trying  to  minimise  the  waveform  error  for  the  entire 
target  sequence.  In  contrast,  a  much  poorer  initial  param¬ 
eter  estimate  converges  to  the  true  parameter  value  for  a 
target  sequence  length  of  87  when  our  modified  minimiza¬ 
tion  method  is  employed. 


dence  into  a  nonlinear  difference  equation  framework.  We 
explained  how  such  sensitivity  could  produce  conditions  ill- 
suited  for  a  proposed  gradient  descent  minimization  of  the 
waveform  error,  and  proposed  an  improved  method  to  over¬ 
come  its  limitations.  The  improved  method  performed  sig¬ 
nificantly  better  than  the  original  when  tested  on  sequences 
generated  from  two  chaotic  dynamical  systems. 
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Figure  4:  Parameter  estimation  performance  for  the 
H»on  system.  Just  as  in  the  one-parameter  case,  the  orig¬ 
inal  minimization  procedure  procedure  falls  immediately 
into  a  local  minimum.  However,  while  the  modified  wave¬ 
form  error  method  reduces  initially  much  larger  parameter 
deviations  better  than  the  original,  it  does  not  converge 
to  the  true  parameters,  due  to  the  singular  nature  of  the 
waveform  error’s  Hessian. 


Figure  6;  Waveform  error  surface  in  the  neighborhood  of 
the  global  minimum.  The  Hessian  appears  to  m  singular. 


